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One of the principle efforts in cosmic microwave background (CMB) research is measurement of 
the parameter /ni that quantifies the departure from Gaussianity in a large class of non-minimal 
inflationary (and other) models. Estimators for /ni are composed of a sum of products of the 
temperatures in three different pixels in the CMB map. Since the number ~ A'^^ix of terms in this sum 
exceeds the number A^'pix of measurements, these ~ A'^^ix terms cannot be statistically independent. 
Therefore, the central-limit theorem does not necessarily apply, and the probability distribution 
function (PDF) for the /ni estimator does not necessarily approach a Gaussian distribution for 
A^'pix 2> 1. Although the variance of the estimators is known, the significance of a measurement 
of /ni depends on knowledge of the full shape of its PDF. Here we use Monte Carlo realizations of 
CMB maps to determine the PDF for two minimum-variance estimators: the standard estimator, 
constructed under the null hypothesis (/ni = 0), and an improved estimator with a smaller variance 
for /ni 7^ 0. While the PDF for the null-hypothesis estimator is very nearly Gaussian when the true 
value of /ni is zero, the PDF becomes significantly non-Gaussian when /ni 7^ 0. In this case we 
find that the PDF for the null-hypothesis estimator /ni is skewed, with a long non-Gaussian tail at 
fni > I /nil and less probability at fni < \fni\ than in the Gaussian case. We provide an analytic fit to 
these PDFs. On the other hand, we find that the PDF for the improved estimator is nearly Gaussian 
for observationally allowed values of /ni. We discuss briefly the implications for trispectrum (and 
other higher-order correlation) estimators. 

PACS numbers: 



I. INTRODUCTION 

The simplest single- field slow-roll inflation models pre- 
dict that primordial perturbations should be nearly 
Gaussian ^ , but with predictably small departures from 
Gaussianity [2] . This is often quantified through the non- 
Gaussianity parameter /„] defined by 3J, 

'i> = + /„i(</)2-((/.2)), (1) 

where $ is the gravitational potential and (j) a Gaussian 
random field. Standard single-field slow-roll infiation pre- 
dicts /ni <C 1 for the primordial field (although nonlin- 
ear evolution of the density field may produce /ni ^ 1 
at the time of recombination; see, e.g., Ref. [4|). How- 
ever, multi-field [5j or curvaton [6] models, or models with 
sharp features [7] or wiggles [S] may produce larger values 
of fn\- Measurement of fni has thus become one of the 
primary goals of cosmic microwave background (CMB) 
and large-scale-structure (LSS) research. Current limits 
from the CMB/LSS are in the ballpark of |/ni| < 100 
HH]- The plot has thickened with a suggestion |11J 
(not universally accepted) that WMAP data prefers (at 
the 2.8(7 level) /„! 0, with a best-fit value fni — 35. The 
Planck satellite [T2] is expected to achieve a sensitivity 
of fni ~ 5. 

In this paper, we address the following question: What 
is the probability distribution function (PDF) P{fni) for 
an estimator /„! that is constructed from a CMB map? 
If the PDF departs from the Gaussian distribution that 
is often assumed, then the 99.7% confidence level (C.L.) 



interval for fni may be different than three times the 
standard deviation for fni- The interpretation of mea- 
surements thus requires knowledge of this PDF. 

The question arises as the theory predicts not only 
the mean value of the estimator fni, but it also makes 
a prediction for the detailed functional form of the PDF 
P(/„i). The consistency of a given measurement of fni 
with a theoretical prediction for fni depends on knowl- 
edge of the shape of P{fni). Thus, for example, we often 
evaluate or forecast the standard error crj^j with which 
a given measurement will recover the true value of 
and then simply assume that the error is Gaussian. If 
so, then with cr/^j = 10, for example, a measurement of 

fni = 30 would represent a 3cr departure from fni — and 
a measurement fni — would represent a 3a departure 
fom fni = 30. However, if the PDF depends on the true 
value fni, and if that distribution is non-Gaussian, then 
it may be that a measurement fni = 30 could be easily 
consistent with a true value fni = 0, while a measure- 
ment fni = could be inconsistent with fni = 30 with 
a confidence greater than "3a." We will see below that 
something like this actually occurs with measurements of 

fni- 

This question is particularly important for measure- 
ments of non-Gaussianity (as opposed, for example, for 
the CMB power spectrum), because fni is a sum over 
products of three temperature measurements (unlike the 
power spectrum, which sums over squares of temperature 
measurements) . Suppose the temperature is measured in 
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A^pix pixels. There are then ~ ^pix terms in the /„! esti- 
mator (after restrictions imposed by statistical isotropy) . 
While these terms may have zero covariance, they are not 
statistically independent; there is no way to construct 
A'^jx statistically independent quantities from iVpix mea- 
surements! The conditions required for the validity of the 
central-limit theorem are therefore not met, and P{fni) 
will not necessarily approach a Gaussian in the A^pix S> 1 
limit. 

The PDF can be obtained from Monte Carlo simula- 
tions, but the simulations are very computationally in- 
tensive (e.g., Ref. [13|). The number of Monte Carlo re- 
alizations is thus usually limited to the number, < 1000, 
required to determine a 99.7% C.L. detection or some- 
times even fewer if it is just the variance that is being 
estimated. Although with only 1000 realizations Fig. 8 in 
Ref. 13] shows hints of a non-Gaussian PDF, simulations 
done up until now do not include enough realizations to 
precisely map the functional form of P(/ni). The num- 
ber of realizations required to map ultimately the 4fT, 
5a, etc. ranges will be prohibitive, especially since the 
simulations will need to be re-run repeatedly to deter- 
mine how the error ranges depend on cosmological pa- 
rameters, instrument-noise properties, scanning strate- 
gies, etc., and they then must be run for multiple theo- 
retical values fni- 

Work along these lines was begun in Ref. [13] , wherein 

it was shown that the variance of the distribution P{fn\) 
may have a strong dependence on the true underlying 
value of /ni- More precisely, they evaluated the vari- 
ance of the estimator designed to have the minimum vari- 
ance under the null hypothesis /„i = (which we refer 
to frequently below as the "null-hypothesis minimum- 
variance" estimator, or NHMV estimator), and showed 
that the variance of this NHMV estimator increases as 
/ni^ increases. They then constructed an alternative es- 
timator f",, which we call the CSZ estimator^, which 
has a PDF with a variance that saturates the Cramer- 
Rao bound up to corrections of order /ni^. Still, as we 
have argued above, the consistency of a hypothesis with 
a measurement requires full knowledge of the PDF of 
whatever estimator is used in the analysis. 

To address these questions, we calculate the PDF for 
an ideal (no-noise) map to understand the irreducible 
PDF introduced simply by cosmic variance under the 
Sachs- Wolfe approximation and on a flat sky. We hope 
that lessons learned about P{fni) in this ideal case may 
help interpret and understand current/forthcoming re- 
sults and assess the validity of full-experiment simula- 
tions. 



We calculate these PDFs by using A-Ionte Carlo real- 
izations of numerous no-noise flat-sky CMB maps. The 
first order of business with a map will be to determine 
whether a given map is consistent or inconsistent with the 
null hypothesis /ni = 0. Therefore, we first calculate the 
PDF that arises if /ni does indeed vanish, for the NHMV 
estimator fni, and we also calculate the PDF that arises 
if the true value of /ni is nonzero. We provide an analytic 
fit for these PDFs in Eq. (21 ). If the evidence from such 



a measurement were to show that f^i is nonzero, then the 
next step would be to apply the CSZ estimator for 
fni 7^ [13] to obtain a more precise value for /„! or to 
test consistency of the data with a specific nonzero value 
of /ni- We therefore follow by calculating the PDF for 
these improved non-null-hypothesis estimators. 

We find that, besides having a variance that increases 
with /„i^, the PDF of the NHMV can have a signifi- 
cantly non-Gaussian shape when f^i ^ with a long 
non-Gaussian tail for /„! > |/ni| and less probability at 

/ni < I /nil than in the Gaussian case. As an example, 
taking = 100 for an experiment which measures mul- 
tipoles out to linax = 3000 (such as Planck) and assuming 
a Gaussian PDF for the NHMV this experiment mea- 
sures 74 < < 148 at the 99.7% C.L.; the actual PDF 
shows that this experiment measures 68 < /ni < 143 at 
the 99.7% C.L. Applying the CSZ estimator to the data 
we find it has a PDF which is well approximated by a 
Gaussian with /^i = 100 ± 12.5 at 99.7% C.L. 

This paper is organized as follows. In Sec. [IT] we con- 
struct the standard minimum-variance estimator /„; un- 
der the null hypothesis — and discuss why the PDF 
for this estimator is not necessarily Gaussian, even in 
the limit of a large number of pixels. In Sec. HI A| we use 
Monte Carlo calculations to evaluate the PDF P{fni) for 
this estimator if the null hypothesis is indeed valid, i.e., 
if /ni is indeed zero. We find that the PDF in this /ni = 
case is well approximated by a Gaussian, for Npi^ ^ 1, 
even though the central-limit theorem does not apply. 
In Sec. |III B[ we calculate the PDF assuming that the 
null hypothesis is not valid, i.e., if fn\ ^ 0. We find the 
PDFs in this case can be highly non-Gaussian, skewed to 
large |/ni|, with long large- /ni non-Gaussian tails and less 
likelihood at /„] < |/ni| relative to the Gaussian distri- 
bution of the same variance. We provide fitting formulas 
for the PDF as a function of the estimator /„!, the true 
value of /ni, and the maximum multipole moment Zmax 
of the map. In Sec. IV we discuss the PDF of the CSZ 



^ We note that the CSZ estimator, which is defined under the 
Sachs- Wolfe hmit, has yet to be generalized so that it can be 
applied to actual data. On the other hand a Bayesian approach, 
discussed in Ref. |15) . allows for an /ni inference that saturates 
the Cramer-Rao bound even in the presence of non-Gaussianity. 



estimator. We show that this estimator is well approx- 
imated by a Gaussian for values of /„; still allowed by 
observations. In Sec. |V]we summarize and discuss some 
possible implications of the work for other bispectra and 
also for the trispectrum and other higher-order statistics. 
An Appendix discusses the computational techniques we 
used in order to perform our Monte Carlo simulations. 
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II. NON-GAUSSIANITY ESTIMATORS 
A. Formalism 

We assume a flat sky to avoid the complications 
(e.g., spherical harmonics, Clebsch-Gordan coefficients, 
Wigner 3j and 6j symbols, etc.) associated with a spheri- 
cal sky, and we further assume the Sachs- Wolfe limit. We 
denote the fractional temperature perturbation at posi- 
tion on a flat sky by T{d) and refer to it hereafter simply 
as the temperature. 

The field T{9) has a power spectrum C; given by 



TrTf)— rtdf . r „Ci, 
where = 47r/sky is the survey area (in steradian) 



(2) 



Tj^ / d"0e-""r(^) ~ — ^e-''-^r(^), (3) 



pix 



is the Fourier transform of T{9), and ^ is a Kro- 

necker delta that sets li = —l^- The power spectrum for 
T{0) is given by 



Ci = 



2'kA 

1^' 



(4) 



where the amplitude, A ~ 10 . The bispectrum 
B{li,l2,h) is defined by 



The Kronecker delta insures that the bispectrum is de- 
fined only for -|- Z2 -l- ?3 = 0; i.e., only for triangles in 
Fourier space. Statistical isotropy then dictates that the 
bispectrum depends only on the magnitudes li, I2, I3 of 
the three sides of this Fourier triangle. 



B. The null- hypothesis minimum- variance 
estimator 

We now review how to construct the minimum- 
variance estimator for /„i under the null hypothesis. This 
is the quantity that one would first determine from the 
data to check for consistency of the measurement with 
the null hypothesis fni = 0. 

From Eq. each triangle h + h + h = gives an 
estimator, 



(/nl) 



123 



Tj- Tj- Tj- 

'■I 1-2 ^3 



(6) 



and under the null hypothesis this has a variance propor- 
tional to 



[nB{hMM)/U 



(7) 



The null-hypothesis minimum- variance estimator is con- 
structed by adding all of these estimators with inverse- 
variance weighting. It is [121 E] 

^_ 2 V- T^TjTjB{hMM)lfni 



and it has inverse variance. 



E 



Il+l2+l3=0 



[B{hMM)/Li 



(9) 



C. Non-gaussianity of the PDF 

If the number of pixels in the CMB map is N^i^, then 
there are also iVpix statistically independent Tp But there 
are a much larger number, oc iVpj^, of triplets T^-^Tj-^Tj-^, 
included in the estimator [cf., Eq. ([s])], and so the num- 
ber of individual "data points" (i.e., triplets) used in 
the minimum- variance estimator scales like iVpi^ 3> iVpix. 
Since the number of terms included in the estimator is 
greater than the number of independently measured data 
points the standard central-limit theorem does not apply. 
Thus, we cannot assume that the PDF of the estimator 
will approach a Gaussian in the iVpix — >■ 00 limit. 

This contrasts with the estimator C"; cx ^ |T^p of the 
power spectrum C/. While the PDF for C; is not neces- 
sarily Gaussian (it has a X2;+i distribution), it is the sum 
of the squares of statistically independent quantities. The 
central-limit theorem therefore applies, and the distribu- 
tion for Ci does indeed approach a Gaussian for large 
The problems we address here for /„! estimators paral- 
lel those discussed in the literature for the quadrupole 
moment C2, as the distribution for quadrupole- moment 
estimators will be highly non-Gaussian and will also de- 
pend on the underlying theory (see, e.g., Ref. pE]). 



III. THE PDF OF /nl FOR THE LOCAL MODEL 

We now restrict our attention to a family of non- 
Gaussian models in which the temperature T{9) has a 
non-Gaussian component; i.e., 

T0) ^ tie) + 3u{m]^ ~ i^mf)). (10) 

where t(ff) is a Gaussian random field with a power spec- 
trum Ci given in Eq. (|4|. To zero-th order in /„i, the 
power spectrum and correlation function for T(ff) are the 
same as those for t{Q\ Note that T(ff) is, strictly speak- 
ing, the temperature fluctuation, so iT{Q')^ = = Tf^p. 
The bispectrum for this model is 

5(^1,^2,^3) =6/„i(G,Ci, +G,G3 -fC/,G3). (11) 
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The temperature Fourier coefficients can be written 



Tj^=tj^+ fnlStj with 



(12) 



Formally, the sum goes from < l^j < oo, but for a finite- 
resolution map, the sum is truncated at some ^max such 
that the number of Fourier modes equals the number of 
data points. 

We now proceed to evaluate P{f ni', fni, ^max), the PDF 
that arises if the true value is /ni for the NHMV estimator 
/ni and for a map with ^max- To do so, we generated large 
numbers of Monte Carlo realizations of maps according 
to Eq. (12 1, for some assumed value of /„!, and then 
applied the estimator in Eq. ([8| to these maps. Each 
map is simulated in harmonic space from /min = 2 up to 
a maximum multipole Zmax- In order to produce a large 
number of realizations we re-expressed the generation of 
maps and implementation of the estimator in terms of 
fast Fourier transforms as discussed in Appendix A. 



A. The PDF of the null hypothesis 
minimum-variance estimator with /ni = 




FIG. 1: Numerical evaluations of -P(/ni; /ni ~ 0,Zmax). The 
left (right) two panels show the PDF for Zmax = 5 and 
'max = 25 for 10® realizations for a scale-invariant power spec- 
trum. In all panels the PDF has been normalized to have a 
unit variance, and the corresponding Gaussian PDF (with 
the same variance) is shown as the red dashed curve. As Zmax 
gets larger, the PDF tends towards a Gaussian. This is not 
guaranteed by the central-limit theorem since the majority of 
the terms that appear in the estimator are not statistically 
independent. 



First we consider the shape of P{fn\', fni — 0, Zmax), 
the PDF for the NHMV estimator in Eq. @ apphed to a 
purely Gaussian (/„! = 0) map. To do this we generated 
10^ Gaussian realizations and applied the estimator in 
Eq. ^ to generate a histogram of values of /„i. From 

this histogram we determined P{fni', fn\ = 0, ^max) out to 
four times the root- variance, as shown in Fig[T] 

First we note that our simulations verify that the vari- 
ance of the distribution for the null case is well approxi- 
mated by the analytic expression ,221 [T7] , 

SAl"^^^ ln(/inax) 

Additionally our simulations show that out to at 
least four times the root- variance, the PDF P(/ni; fni — 
0, Zmax) is well approximated by a Gaussian for Zmax ^ 25, 
even though the conditions for the central-limit theo- 
rem to apply are not satisfied. Therefore, a measure- 
ment of fni that differed from at more than three times 
the root-variance would indeed constitute a '99.7% con- 
fidence level' inconsistency with the fni ~ hypothesis. 



B. The PDF of the null hypothesis 
minimum-variance estimator with /ni 7^ 

We now consider the form of P(/ni; /ni, 'max) when 
/ni 7^ 0, the PDF for the null- hypothesis minimum- 
variance estimator if the null hypothesis is in fact not 



valid. In this case, the non-Gaussian statistics of the Tjs 
impart some non-Gaussianity to the fni PDF. 

In Fig. [2] we show P{f ni] fni, 'max) calculated using 10^ 
realizations with fni = 1500 and 'max = 25. Clearly the 
PDF in this case is highly non-Gaussian. 

Non-Gaussianity of P(/ni; fni, 'max) for a central value 
/ni 7^ may be significant for the interpretation of data. 
Suppose, for example, that a CMB measurement returns 
fni ~ with a root-variance u/^j = 40. If the PDF 
was assumed to be Gaussian the measurement fni — 
would rule out fni = 100 at the 2.5a level, but given the 
asymmetric PDF of Fig. [2] it may rule out f^i = 100 at a 
much higher significance. 

In order to better understand the origin of the non- 
Gaussian PDF, it is useful to expand the minimum- 
variance estimator in Eq. ([8| to linear order in fni jl4] : 

U1-E0 + fnlEl + ■■■ , (14) 

where 

ll+h+l3=0 

Since Eq ~ and Ei ~ t^, it is clear that (Eq) = 
and (EqEi) = 0, and the normalization guarantees that 
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FIG. 2: The PDF P(/ni) when fni = 1500 using the estimator 
in Eq. ||8| with /max ~ 25. The upper (lower) panel shows 
the PDF on a linear (log) scale. We can see that the PDF 
is significantly non-Gaussian with an exponential drop-off to 
the left of mean and a power-law to the right. We provide a 
fitting formula for P{fnV, fni, /max) in the text. 



(El) = 1. Furthermore, since we have already estab- 
lished that P{fn\) approaches a Gaussian in the large 
^max limit if /ni = 0, we know that, to leading order, 
the non-Gaussian shape of P(/ni; /ni, ^max) f^or /„! 7^ is 
being generated by Ei. 

Some of the statistics associated with Ei have al- 
ready been explored in Ref. [H]. There it is noted 

that the variance of /ni is dominated by Ei in the 
high S/N limit leading to a slower scaling of the S/N 
than the ^^^^ ln~^(^max) scaling expected if the estima- 
tor saturated the Cramer-Rao bound ^4j. We explored 
the same limit using our Monte Carlo realizations, as 
shown in Fig. [3] and find the same qualitative trend but 
with a different dependence on Zmax- Ref. |T3] found 
((Ai?i)^) oc In^^(^niax) whereas our simulations show 
{{AEi)^) cx ln^'^(^i„ax)- We have checked the scaling 
found with our simulations by computing the variance 
analytically, as we further discuss in Appendix B. Fig. [3] 
shows the agreement between our analytic calculation 
(solid curve) and simulations (data points). 

Our simulations allow us to generate the full PDF for 
El, not just the variance. Fig. |4] shows this PDF for 



FIG. 3: The dependence of {(AEi)^) on /max- The points 
correspond to the results of our Monte Carlo simulations for 
1000 realizations at different values of /max. The solid curve 
shows the analytic calculation of the variance presented in 
Appendix B which is well- fit by the function ((AiJi)^) = 

[14.0(/max)°-^'']/[ln5-^(/max)] ~ 4.5 lu'^ (/max ) . 
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0.0001 



-4 -2 2 4 

(El - {Ei))/aE, 

FIG. 4: The PDF of Ei calculated using lO'* reafizations. 
The thin solid curves correspond to P{Ei) with /max = 25 
(green), /max = 50 (purple), /max = 100 (black). Since the 
functional form of the PDF for each choice of /max is nearly 
identical, we conclude that P{Ei) approaches an asymptotic 
functional form in the large-/max limit. The thick red dashed 
curve corresponds to a fit to P(i?i), accurate to ~10% out to 
3 times the root variance, using the fitting formula in Eq. 1 17 1 
with parameter values Xp — —0.22, a — 0.80, and c = 0.91. 
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various choices of /,nax (thin sohd hnes). An important 
conclusion from Fig. |4] is that the shape of the PDF ap- 
proaches a universal form in the Zmax ^ 1 limit. We 
provide a fit to the PDF (thick red dashed line), accu- 
rate to ^10% (40%) out to three (four) times the root 
variance, using the fitting formula 



log[F(x)] = iV - 



(17) 



X X jj 



where N = y/2/TTa + cexplc^ /a'^]Ki{c^ /a'^) and Ki{x) 
is a modified Bessel function of the first kind, c quantifies 
the non-Gaussianity of the distribution (and approaches 
a Gaussian in the c — )■ cx) limit) and Xp is the value of 
{El — {Ei))/aE^ at the peak of the distribution. The 
red curve in Fig. |4] shows Eq. (171 with parameter values 
Xp = -0.22, a = 0.80, and c = 0.91. 

We are now in a position to write down a semi-analytic 
expression for -P(/ni; /ni, ^max), accurate to ~10% (40%) 
out to three (four) times the root variance, as a function 
of /ni and Zmax- Letting ctq and ai denote the standard 



deviations of the distributions for Eq and Ei respectively 
we have 



1 



8^^max lll(Zmax) 



9/nl 



21n-'(;„ 



(18) 
(19) 



A good approximation to the PDF of /„i is provided by 
the convolution of the PDF of Eq and /ni-Bi : 



P{fnh fn\, I 



(20) 



9v 27rcroO'i 

X / Go{fn\- X)F{[X - fn\]/(Tl)dx 



where Go (a;) is a Gaussian with zero mean and standard 
deviation tJo and F{[x — f^il/ai) is given by Eq. (17 1 with 
Xp = -0.22, a = 0.80, and c = 0.91. 

To obtain an analytic expression for the PDF we can 
approximate the convolution in Eq. (21 1 to write 



n 



afa^ 1 



9 
erf 



exp 



ca'l 



x/2 



G\a 1 + erf 



exp 



erg +criCr^(A: + 0-i) 
X^ 



(21) 



where AT = /„! + XpOx - fn\- 

Another useful way of quantifying the non-Gaussian 
shape of ^'(/ni; /ni, ^max) IS to mcasurc its skewness, 
((A/ni)'^), as a function of f^i and ?max- We show this in 
Fig. [5] for /ni = 100. An analytic fit to the skewness is 
given by 



((A/„i 



Jnl 



/ill 

100 

1 



(22) 



-0.26 



,l + 3.7exp [-(?„ax- 5.1)/740] 
with the variance of the distribution, cr^ ^, given by 

1 



8^^niax In(Zinax) 



1 



max 
In^(^max) 



(23) 



Finally, we note that the shape of i^(/ni; /nh ^max) de- 
parts significantly from a Gaussian when ao — cri. This 
occurs when 



6/n 



(24) 



Therefore, for the Planck satellite (i.e., ^^ax — 3000) the 
non-Gaussian features of P{fn\', fn\, 'max) for the NHMV 
estimator are significant if /„! > 0(10). Thus, given that 
Planck is expected to measure /„; with a variance a ~ 5, 
these PDFs may need to be taken into account to assign 
a precise confidence region with Planck data. 



IV. THE PDF OF AN IMPROVED ESTIMATOR 
WHEN /ni ^ 

As we saw in the previous Section the standard 
(null-hypothesis) minimum-variance estimator /„! is con- 
structed under the null hypothesis, so its variance is 
strictly minimized only when applied to maps with 
/ni = JT^. In particular, the variance of /ni is given 
in Eq. g so that when 36 AJ J l^^J \n^{U,^) > 1, 
the variance scales as the In ^{l-max), as opposed to 
l~^x lii~^('max)- This indicates that when /ni 7^ there 
may be other estimators with smaller variances. 

For a flat-sky and under the Sachs- Wolfe approxi- 
mation Ref. |14| introduced an improved estimator for 
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1000 2000 3000 4000 

^max 

FIG. 5: The skewness, {(A/ni)^), as a fraction of the variance 
of P(/ni; /ni, Zmax) as a function of Zniax for /nl = 100. We 



provide an analytic fitting formula in Eqs. (221 and (231 as a 
function of f^i and /max- 



/nl 7^ which has a variance that continues to decrease 
as l/[Z^ax lii(^max)] in the high signal-to-noise limit. 
To achieve this scaling they introduced a realization- 
dependent normalization, 



Jnl 



E 



where 



h+l2+h=0 



xr 



(25) 



(26) 



By construction {M) — 1. They then define a new esti- 
mator constructed under the non-imll hypothesis: 



f n — 
Jnl — 



nl 



(27) 



To explore the properties of the PDF of we expand 
the normalization as A/" ~ A/q + /niA/i + • ■ ■ and write 



f n 
Jnl 



Eo 



f/nl- 
./nl^l 



EiAfo - EnNi 



A/? 



(28) 
(29) 



In order to determine the shape of P{fni)i we com- 
puted P{£o) and P{£i) for various values of ^,nax- We 
found, as in the case, that these PDFs approach 
asymptotic shapes in the Zmax S> 1 hmit. We show 
these PDFs in Fig. [g] determined by 10^ realizations for 
^max = 25. It is clear that P{£o) is very well approxi- 
mated by a Gaussian, whereas P{£i) has significant non- 
Gaussian wings. As in the P(/ni) case, this implies that 




-2 s 

(fo - (fo))/%, 



-2 2 



FIG. 6: The PDF P{£o) (left) and P{Si) (right) for Z^ax = 25 
determined with 10^ non-Gaussian realizations. The top pan- 
els show the PDF on a linear scale; the bottom panels show 
the PDF on a log scale. We have confirmed that the shape 
of the PDF is unchanged for larger values of Zmax. The PDF 
of £o (left) is well approximated by a Gaussian. However, 
the PDF of the first-order correction £i (right) has significant 
non-Gaussian wings. This implies that the full PDF of is 
also non-Gaussian, even if the true value of fni matches that 
assumed in the construction of the CSZ estimator. Quantita- 
tively, however, the level of non-Gaussianty will be small for 



Planck, as the variance of £i is ((Afi) 



91n^(Zn,ax)/(4ax). 



the level of non-Gaussianity in Pif^^^i) is significant only 
when the ratio /„i^ {{^£if) I {{^£of) > 1- Our simu- 
lations show 



((A^o)^) 
((Afi)^) 



1 



8^^max 



In (^max 



ln(Zn 

) 



c)' 



(30) 
(31) 



so that the PDF will be significantly non-Gaussian when 



/nlAl/2 > l 



81n(/n 



1/2 



(32) 



Therefore, for Planck (with /„,ax = 3000) P(/« ; /„i, U^^) 
will be significantly non-Gaussian only if /„i > 0(1000). 
Since this has already been ruled out by observations [21 
[TU] . we conclude that -P(/^';; /ni, ^max) will be effectively 
Gaussian. 



V. DISCUSSION 

Here we have argued that the PDF for non-Gaussianity 
estimators cannot be assumed to be Gaussian, since the 
number of triplets used to construct these estimators may 
greatly exceed the number A^pix of measurements. The 
99.7% confidence- level interval cannot safely be assumed 
to be 3 times the 66.5% confidence- level interval. We 
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found, however, that the standard minimum- variance es- 
timator /ni constructed under the null hypothesis is well- 
approximated by a Gaussian distribution in the Zmax S> 1 
limit if the null hypothesis is correct (i.e., when applied 
to purely Gaussian maps). 

We then calculated the same PDF -P(/ni) under the 
hypothesis that the true value of fni is non-zero. We 
find that the PDF is non-Gaussian in this case, skewed 
to large /„! if /ni > and vice versa for < 0. The 
PDF for small positive or for negative fn\ is significantly 
smaller for /„i > than the Gaussian PDF with the same 
variance. Thus, for example, if the NHMV estimator 
gives /ni > 0, it may actually rule out /ni = with 
a smaller statistical significance than would be inferred 
assuming a Gaussian distribution of the same variance. 
For Planck (with Zmax — 3000) we find that the non- 
Gaussian shape of P(/ni) is significant if /„i > 0(10). 
Thus, the non-Gaussian shape of the PDF may need to 
be taken into account, even in case of a null result, to 
assign a precise 99.7% confidence-level upper (or lower, 
for /ni < 0) limit to We also provide, in Eq. (21 1, an 
analytic fit to these PDFs. 

The non-Gaussian shape of P{fni) when /„! 7^ is 
accompanied by a variance that decreases only logarith- 
mically with increasing /max- Because of this, Ref. [13] 
constructed an improved estimator under the /„! ^ hy- 
pothesis with a variance that saturates the Cramer-Rao 
bound and continues to decrease as l/[lf-^^^\og{lniax)]- 
We found that for observationally allowed values of /ni 
this improved estimator has a PDF that is well approx- 
imated by a Gaussian shape. However, this estimator 
has only been defined under the Sachs- Wolfe limit and it 
is not immediately clear how it should be generalized to 
be applied to actual data. An alternative, Bayesian, ap- 
proach to measuring f^i which also saturates the Gramer- 
Rao bound in the presence of /ni 7^ is presented in 
Ref. [H]. 

We have restricted our attention to the bispectrum in 
the local model, but the PDF must be similarly deter- 
mined for the non-Gaussianity parameter for bispectra 
with other shape dependences; e.g., the equilateral model 
[iniHn' or that which arises with self-ordering scalar fields 
[5T] . It should also be interesting to explore the PDF for 
maximum-likelihood, rather than quadratic, estimators 
(see, e.g., Ref. [H]). Ultimately, a variety of experimen- 
tal effects and more precise power spectra and bispectra, 
rather than the Sachs- Wolfe-limit quantities used here, 
will need to be included in interpreting the results of re- 
alistic experiments. 

There is also interest in using higher-order correlation 
functions to measure f^i from CMB maps. Our argu- 
ments should apply also to these higher-order correlation 
functions, like the trispectrum, etc. For example, the 
estimator for the amplitude of the n-point correlation 
function (e.g., n ~ 3 for the bispectrum, n — 4 for the 



trispectrum, etc.), will be constructed from . .p^^ 



N. 



(n-l) 



tions scales even more rapidly with A'pix than that for 
the bispectrum. Thus, although the signal-to-noise scales 
more rapidly with iVpix for these higher-order correlation 
functions than that for the bispectrum [T71 [221 [53] , con- 
cerns about the PDF for these estimators should be even 
more serious than for the bispectrum. It will thus be 
necessary to understand the PDF for these higher-order 
estimators to confidently forecast the statistical signfi- 
cance of measurements 1241. 
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Appendix A: Computing non-Gaussianity estimators 
using FFTs 

We are interested in using Monte Garlo simulations to 
determine the shape of the PDF of fni as a function of 
the fiducial choice of /ni and the number A'pix of pixels 
measured in a given observation. Applying the estimator 
in Eq. ^ to the local-model bispectrum [Eq. ( 11 )] it can 
be rewritten 



/nl 



Jnl 



E 



|il+i2+/3|=0 



Tj- Tj- Tj- 

tl t2 ^3 



(Al) 



The estimator in Eq. (Al) takes N^^^ operations to eval- 
uate. Since current GMB observations have A^pix ~ 10^ 
this estimator would take a prohibitively long time to 
evaluate for a significant number of realizations, espe- 
cially since we are interested in probing the shape of the 
PDF far into the tail of the distribution (~ 3 — 4cr). 

As discussed at length in Ref. [IS] this is even more of 
a problem when measuring non-Gaussianity on the full 
sky where the number of operations scales as N^l^ . In 
order to make the problem tractable Ref. [53] rewrites 
/nl in terms of real-space quantities reducing the number 

3 /2 

of operations to N^l^ . 

We can do the same for /„! in the flat-sky approxima- 
tion. Noting that 



combinations of n pixels, and this number of combina- 



n 



pie-ih+h+h) 



(A2) 
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and writing 



[Eq. ([16])] defined by 



r 



/ni can be written 



/nl 



Jnl 



IT 



^2(e»)B(6i). 



(A3) 
(A4) 



(A5) 



Next, in order to compute the integral in Eq. ( |A5[ ) we 
use the Nyquist sampUng theorem and the fact that both 
A{9) and B{d) have finite Fourier spectra (truncated at 
a maximum frequency Zmax)- This allows us to rewrite 
the integral as a discrete sum 



/nl 



rial ^ 



Jnl 



X B\2ti 



i=i j=i 

i 



2n- 



N 



1 9^^'-l 

N 



N ' 



1 „ i - 1 

N 



(A6) 



where iV = 2{2l„ 



!)■ 



Since Eqs. (A3) and (A4) are discrete inverse Fourier 



transforms we can use a fast Fourier transform (EFT) 
algorithm so that the number of operations scale as 
A^pix In(A^pix). 

We can use the same computational trick when evalu- 
ating the non-Gaussian contribution for each realization 
by also employing a forward EFT in order to compute 
the convolution in Eq. (12 1. 



Appendix B: Analytic calculation of ((A_Bi)'^) 

In order to verify that our simulations are correct we 
performed an analytic calculation of the variance of Ei 



li h 13 



B{hMM)- (Bl) 



li+l-2+h=0 



A straightforward but tedious calculation shows that the 
variance is given by 



where 



{{^Eif) = 9a^_J^i + 8^2 + A3 + 4^4), (B2) 



Aa 



B{1) B{k) 



Sr 



{/>,{fc} 

B{l)B{k) 

{i},{k} 

{1} '1 ' ' = l 

^ B{l)B{k) ^ 

{r),{k} ' ' ' 



(B3) 
(B4) 

(B5) 
(B6) 



where {/} indicates the sum is over \li + h + hi =0 and 
B{1) = B{li, I2, h). Computing these terms as a function 
of /max we find that the variance is well-fit by the function 



70.433 



(B7) 



In Fig. [3] we show how that this analytic calculation of 
the ((Ai5i)^) is reproduced by the results of the Monte 
Carlo simulations. 
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